#preparing data for Stu

#1 - all sequences and snps in clustal format

clustalOutput <- function(seq.name,dir.embl,snp.df,dir.out,name.out=NULL){
	#create clustal alignments (as a series of fasta sequences)
	#requires cDNA name, directory of embl files, df containing snps (grouse.df3)
	#an output directory and, optionally, an output file name
	#much of code hacked from coding_snps3.R
	
	snp.df <- snp.df[snp.df$cDNA==seq.name,] #note assumes column cDNA
	if(is.null(snp.df)|nrow(snp.df)==0) return('No SNPs found for ',seq.name,'... returning\n')	
	embl.file.name<-paste(dir.embl,"/",seq.name,"_exon.embl",sep="")

	if(!file.exists(embl.file.name)){
		error(embl.file.name, " does not exist... returning\n")
		}
	
	#read sequence
	ref.seq <- readEmblSeq(embl.file.name)

	#find masks
		mask.lines <- grep("/note=\"mask\"",tmp.features) -1
	if(length(mask.lines)>0){
		grep.mask <- tmp.features[mask.lines]
		grep.mask <- sub('^[ ]*misc_feature','',grep.mask)
		grep.mask <- gsub('[ ]*','',grep.mask)
		
		mask.list <- strsplit(grep.mask,"\\.\\.")
		mask.list <- lapply(mask.list,as.numeric)
		mask.list <- lapply(X=mask.list,FUN=function(X){
			if(length(X)==1){
				return(c(X,X))
				}else{
				return(X)
				}
			})
		mask.list <- matrix(unlist(mask.list),ncol=2,byrow=T)
		colnames(mask.list) <- c("start","end")
		mask.list <- as.data.frame(mask.list) 
		#output starts and ends of mask as a dataframe
		tmp.loc <- numeric(0)
		for(i in 1:nrow(mask.list)) tmp.loc <- c(tmp.loc,seq(mask.list[i,1],mask.list[i,2]))
		for(i in 1:length(tmp.loc)){
			ref.seq <- paste(substr(ref.seq,1,(tmp.loc[i]-1)),"-",substr(ref.seq,(tmp.loc[i]+1),nchar(ref.seq)),sep="") #using "-" as mask for clustal
			#will break if masked nucl at end of sequence
			}

		}#end of mask find
	
	#mutate sequence
	mut.seq <- character(length=nrow(snp.df))
	for(snp in 1:nrow(snp.df)){
		tmp.id <- strsplit(snp.df[snp,"id"],"_")[[1]]
		new.seq <- tmp.id[length(tmp.id)] #picks off last part of id
		mut.seq[snp] <- paste(
			paste(rep("-",snp.df[snp,"start.pos"]-1),collapse=""),
			new.seq,
			paste(rep("-",nchar(ref.seq)-(snp.df[snp,"start.pos"]-1)-nchar(new.seq)),collapse=""),
			sep="")
		}
	if(is.null(name.out)) name.out <- paste(seq.name,"clustal.fasta",sep="_")
	name.out <- paste(dir.out,name.out,sep="/")
	cat(paste(">",seq.name,"reference",sep="_"),
		paste(ref.seq,"\n",sep=""),
		file=name.out,append=FALSE,sep="\n")
	for(snp in 1:length(mut.seq)){
		cat(paste(">",seq.name,snp.df[snp,"start.pos"],sep="_"),
			paste(mut.seq[snp],"\n",sep=""),
			file=name.out,append=TRUE,sep="\n")

		}
	}

#debug(clustalOutput)
#clustalOutput(seq.name="cDNA_271-1",dir.embl="../embl_exon3",snp.df=grouse.df3,dir.out="../test_files")
dir.create('/Paterson/Datafiles/grouse/clustal_snps')
for(cdna in unique(grouse.df3$cDNA)){
	clustalOutput(seq.name=cdna,dir.embl="/Paterson/Datafiles/grouse/embl_exon3",snp.df=grouse.df3,dir.out="/Paterson/Datafiles/grouse/clustal_snps")
	}

#2 - export allele frequencies etc into excel
write.table(cbind(grouse.df3[,c("cDNA","id","start.pos","pop.1.freq","pop.2.freq", "pop.3.freq", "pop.1.tsv.depth", "pop.2.tsv.depth","pop.3.tsv.depth", "Fst", "pop.1.count", "pop.2.count", "pop.3.count","mean.het")],BBlikTest.df),quote=FALSE,row.names=FALSE,sep="\t",file="/Paterson/Datafiles/grouse/Stu_output/allele_freqs.txt")
write.table(BBlikTest.pooled.df,quote=FALSE,row.names=FALSE,sep="\t",file="/Paterson/Datafiles/grouse/Stu_output/fst_pooled.txt")
write.table(grouse.chick.MK[,c("ensembl.id",  "grouse.id", "ensembl_transcript_id", "Dn", "Ds","Pn", "Ps", "fisher.pv", "fisher.or", "fisher.lower",  "fisher.upper","transitions","transversions","dN", "dN.se", "dS", "dS.se", "kappa","omega", "Ls", "Ln", "t","Fst.mean")],quote=FALSE,row.names=FALSE,sep="\t",file="/Paterson/Datafiles/grouse/Stu_output/divergence_summary.txt")

save.image('wksp_2011_01_11.rdata')

write.table(grouse.TjD.df,file="/Paterson/Datafiles/grouse/Stu_output/TajimaD_summary.txt",,quote=FALSE,row.names=FALSE,sep="\t")